clear;
% 导入CSV数据
filename = '5kA0.csv';
data = readmatrix(filename);
t = data(:, 1);  % 时间列
y = data(:, 3);  % 电流列

% 数据去重：同一时间点可能存在多个电压值，取中值
[uniqueTime, ~, idx] = unique(t);
yClean = accumarray(idx, y, [], @median); % 使用中值减少异常值影响
tClean = uniqueTime; 


% 重新排序数据（若原始时间未排序）
[~, sortedIdx] = sort(tClean);
tSorted = tClean(sortedIdx);
ySorted = yClean(sortedIdx);

figure;
plot(tSorted,ySorted);
% 设置坐标轴范围
axis([-0.001 0.01 -4000 30000]);
grid on;
% 计算均匀时间步长

meanDeltaTime = mean(diff(tSorted)); % 计算平均时间间隔
meanDeltaTime1 = diff(tSorted); % 计算平均时间间隔




% % 创建均匀时间轴（可能根据数据量调整点数，如1e5）
% numPoints = 100000; % 总点数
% tUniform = linspace(min(tSorted), max(tSorted), numPoints)';
% dtUniform = tUniform(2) - tUniform(1); % 实际均匀步长
% 
% % 插值到均匀时间轴
% yUniform = interp1(tSorted, ySorted, tUniform, 'linear', 'extrap');


% 设置滤波器参数（示波器为50MHz采样则dt~20ns, 即fs=1/dtUniform）
fs = 1 / 1e-8;    % 采样率
% 奈奎斯特频率
Fn = fs / 2;

%  滤波器设计：5阶低通 Butterworth 滤波器，截止频率设为有效频率的 0.1 倍 ==
order = 5;              % 滤波器阶数
Fc = 0.2 * Fn;          % 截止频率（例子中为有效频率的10%）
cutoffNorm = Fc / (fs/2); % 归一化截止频率 

% 设计巴特沃斯滤波器
[b, a] = butter(5, cutoffNorm, 'low');

% 零相位滤波（filtfilt消除相位畸变）
filteredY = filtfilt(b, a, ySorted);




% 检测过零事件
zeroCrossings = filteredY(1:end-1) .* filteredY(2:end) <= 0; % 找到符号变化区域
zcIndices = find(zeroCrossings)';

% 计算每个过零点的准确时间
zcActualTimes = [];

for zci = zcIndices
    t1 = tSorted(zci);
    t2 = tSorted(zci + 1);
    v1 = filteredY(zci);
    v2 = filteredY(zci + 1);
    if v1 < 0 && v2 > 0 || v1 > 0 && v2 < 0
        % 线性插值确定过零时间  
        tzc = t1 - v1 * (t2 - t1) / (v2 - v1);
        zcActualTimes = [zcActualTimes; tzc];
    end
end


% figure; 
% plot(tUniform, filteredY);
% hold on;
% plot(zcActualTimes, zeros(size(zcActualTimes)), 'rx'); % 过零点标记
% title('滤波后导数及过零标记');
% hold off;
% grid on;

window = 1e-6; % 过零点前后各2us
fitOrder = 2; % 可调整的多项式阶数（尝试用3次曲线）
slopes = [];
tzc_final = [];
for tzc = zcActualTimes'
    if tzc>0.003 && tzc<0.006
        % 时间窗口范围
        windowStart = tzc - window;
        windowEnd   = tzc + window;
        
        % 筛选出时间窗口对应的索引
        windowMask = (tSorted >= windowStart) & (tSorted <= windowEnd);
        t_window   = tSorted(windowMask);
        v_window   = filteredY(windowMask);
    
        % 避免窗内数据点太少（例如少于阶数+1）触发错误
        if length(t_window) < fitOrder + 1
            warning('窗口内的数据点不足以进行%u阶多项式拟合，请调整采样率或窗口范围.', fitOrder);
            continue;
        end
    
        % 移动时间基到窗口中心，提高拟合精度
        % tNorm = t_window - tzc;
        % [p_fit,~,mu] = polyfit(t_window, v_window, fitOrder);  缩放之后求导错误
        p_fit= polyfit(t_window, v_window, fitOrder);
        % 拟合导数：dp = p(1)*fitOrder + ...  或者使用MATLAB的导数函数
        dpols = polyder(p_fit); % 原始阶数的导数多项式
        slope_zc = polyval(dpols,tzc)/1e6; % 过零点导数
        if slope_zc <-100
            tzc_final = [tzc_final;tzc];
            slopes = [slopes; slope_zc];
        end
    end
end



figure; 

plot(tSorted*1e3, filteredY/1e3);
hold on;
plot(tzc_final*1e3, zeros(size(tzc_final)), 'gx'); % 过零点标记
%设置全局字体，放在plot之后
set(gca,'FontName','Times New Roman','FontSize',18);  
%单独设置title字体，需放在set（gca）之后
title('滤波后导数及过零标记','FontName','黑体','FontSize',16); 
hold off;
grid on;
%设置坐标轴label
xlabel('time (ms)');
ylabel('Current /kA');

% 设置坐标轴范围
axis([-2 10 -10 10]);

fprintf("%s 过零点%.4f ms的di/dt为 %.2f A/us. \n", filename,tzc_final(1)*1000,slopes(1));

